Generic properties of a quasi-one dimensional classical Wigner crystal 
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We studied the structural, dynamical properties and melting of a quasi-one-dimensional system 
of charged particles, interacting through a screened Coulomb potential. The ground state energy 
was calculated and, depending on the density and the screening length, the system crystallizes in 
a number of chains. As a function of the density (or the confining potential), the ground state 
&X), configurations and the structural transitions between them were analyzed both by analytical and 

^ ■ Monte Carlo calculations. The system exhibits a rich phase diagram at zero temperature with 

' continuous and discontinuous structural transitions. We calculated the normal modes of the Wigner 

crystal and the magneto-phonons when an external constant magnetic field B is applied. At finite 
temperature the melting of the system was studied via Monte Carlo simulations using the modified 
Lindemann criterion (MLC). The melting temperature as a function of the density was obtained 
for different screening parameters. Reentrant melting as a function of the density was found as 
well as evidence of directional dependent melting. The single chain regime exhibits anomalous 
' melting temperatures according to the MLC and as a check we study the pair correlation function 

(~ i ' at different densities and different temperatures, which allowed us to formulate a different criterion. 

, Possible connection with recent theoretical and experimental results are discussed and experiments 

^ ■ are proposed. 



I. INTRODUCTION 



Recently there has been a great deal of interest in niesoscopic systems consisting of interacting particles in low 
'■^ . dimensions or confined geometries. A class of quantum anisotropic systems exhibiting "stripe" behavior appears in 
^ ' the quantum Hall effect ]3| , in oxide manganites and high-Tc superconductors Q where electronic strong correlations 
O I are responsible for the formation of these inhomogeneous phases. Another class of confined quasi one-dimensional 
(QID) geometries appears in many diverse fields of research and some typical and important examples from the 
experimental point of view are : electrons on liquid Helium i^] , microfluidic devices ||5| , colloidal suspensions 
and confined dusty plasma 0. 

^ ' A major phenomenon which is expected to occur in charged particles interacting via a Coulomb or screened Coulomb 
potential is Wigner crystallization (WC) !3| at low enough temperatures and densities when the potential energy 
overwhelms the kinetic energy. Indeed, evidence of such a type of transition was found very recently Q in experiments 
QQ ' on electrons on the surface of liquid Helium where the electrons were confined by metallic gates and exhibited 
' dynamical ordering in the form of filaments. This particular experiment posed many interesting questions regarding 
' the nature of the transition to WC, its density dependence and the melting. Furthermore, the considered system 
1 has beeriproposed as a possible step towards the realization of a quantum computer with electrons floating on liquid 
Hehum p. 

In this paper, as a first step towards the understanding of the behavior of these systems, we start with a two- 
dimensional system consisting of an infinite number of charged particles and we impose a parabolic confining potential 
^ in one direction. The particles interact with a Yukawa-type potential where the screening length is an external 
^ parameter. Physically, it can be adjusted e.g. by the gate voltage that confines the electrons. The combination of 
Q \ the interaction among particles and the external potential leads to a rich structural phase diagram as a function of 
O ■ the screening length A and the density n of the system. The structural units (at temperature T = OK) are parallel 
^ ' chains of particles the number of which depend on the values of A and n. The transition from one configuration to 
• i-H 1 the other can be obtained via a first or a second order transition. 

■ Before proceeding further, we should comment on the possibility of two-dimensional crystalline order. According 
^ , to the Mermin- Wagner theorem [lol | there is no true long-range crystalline order in two dimensions. However this 
- - ' theorem is only strictly valid when the potential falls off faster than 1/r and in the thermodynamic limit. When the 
same arguments of the theorem are applied to a large but finite system, no inconsistencies arise from the assumption 
of crystalline order. Thus any system that can be studied in laboratory or in computer simulations can exhibit 
crystalline order pj|. On the other hand short-range order is expected to form even in the thermodynamic limit. 

In a related work [T^ which discussed the temperature equilibration of a one-dimensional Coulomb chain, two 
different equilibration temperatures were assigned {T± and T| | ) , reflecting the different behavior of the modes due to 
the strong confinement. 

The WC in strictly one dimension and in the quantum regime was first studied by Schultz |13l] . He found that for 
arbitrarily weak Coulomb interaction the density correlations at wave vector Akp decay extremely slowly (the most 
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slowly decay term is oc exp (— cVmx) . 

Other remarkable work on the quantum transport and pinning in the presence of weak disorder, where it was shown 
that quantum fluctuations soften the pinning barrier and charge transfer occurs due to thermally assisted tunneling, 
is described in ^T^. 

In addition to the structural properties, it is instructive to study the normal modes of these kind of anisotropic 
systems. There are optical and acoustical branches and their number is equal to the number of chains. The acoustical 
modes correspond to motion along the unconfined direction and the optical ones to motion along the confined direction. 
There is softening of an optical phonon at those values of the density for which we have a continuous structural 
transition. We also study the collective excitations in the presence of a constant magnetic field perpendicular to the 
plane of the system. These modes (magnetophonons), can be directly detected experimentally [m.[l7l|. 

Another important aspect of the problem is the melting as the temperature is raised. The mechanisms of melting 
is of great scientific and technological importance. In infinite two-dimensional (2D) systems theory based on 
unbinding of defects, predicts a two-stage melting where the two stages are continuous. Recent theoretical studies of 
melting of colloidal crystals in the presence of a one-dimensional periodic potential [T^j l revealed a number of novel 
phases and the possibility of reentrant melting. These results depend on the commensurability ratio p — a/dext where 
a is the spacing between the Bragg planes of the 2D system and dext the period of the external periodic potential. 
This kind of system was realized experimentally in 2D colloids in the presence of two interfering laser beams. The 
present work is complementary to the work of Frey, Nelson and Radzihovsky in the sense that a single confining 
potential is considered here, which is not repeated in space. Therefore it can be viewed as a study of a focused 
portion of the infinite 2D system, where we pay attention to only one potential trough neglecting the interaction with 
the rest. With respect to the melting, we found the following remarkable results: (i) a phase diagram which exhibits 
reentrant melting behavior as a function of the density where the different configurations are explored, (ii) a regime 
of frustration exists close to the structural transitions, and (iii) there is evidence that the system first melts in the 
unconfined direction and subseque ntly in the direction where it is confined exhibiting a regime similar to the locked 
floating solid regime found in Ref. jl^ • 

The paper is organized as follows. In Sec. II, we present the model and the methods used. In Sec. Ill we study 
the zero temperature phase diagram and properties of the structural transitions. Sec. IV is devoted to the study of 
the normal modes of the system and in the presence of, or without, an external magnetic field B. In Sec. V we study 
the melting and analyze furthermore in some details the problem of the single chain melting. 

Finally we discuss the connections with recent experimental results and suggest experiments where this behavior 
can be observed in Sec. VI. 



II. MODEL AND METHODS 



The system is modelled by an infinite number of classical charged particles with identical charge q, moving in a 
plane with coordinates r = {x,y). The particles interact through a Yukawa potential and an additional parabolic 
potential confines the particle motion in the y direction. The Hamiltonian of the system is given by : 
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where m is the mass of each particle, e is the dielectric constant of the medium particles are moving in, ujq measures 
the strength of the confining potential. The Hamiltonian can be rewritten in a dimensionless form, introducing the 
quantities vq = {2q/meujQy^^ as unit of length and Eq = (jnuJ^q^ j^e^Yl'^ as unit of energy. Then it takes the form 

where H' — H/Eq, k — tq/X and V'' = ~r^/ro. This transformation is particularly interesting because now the 
Hamiltonian no longer depends on the specifics of the system and becomes only a function of the density and the 
dimensionless inverse screening length. The quantities introduced allow us to define a dimensionless temperature 
r = T/To with To = (mwggV2e2)i/3fc-i. 

For the calculations of the ground state energy we used a combination of analytical calculations and Monte Carlo 
simulations with the standard Metropolis algorithm. This recursive algorithm consists in displacing randomly one 
particle and accepting the new configuration if its energy is lower than the previous one; if the new configuration 
has a larger energy the displacements are accepted with probability d < exp{—AE/T), where (5 is a random number 
between and 1 and AE is the increment in the energy. We have allowed the system to approach its equilibrium 
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state at some temperature T, after executing 10^ 10^ Monte Carlo steps. We have used the technique of simulated 
annealing to reach the T = equilibrium configuration: first the system has been heated up and then cooled down 
to a very low temperature. In the simulations typically 300 particles were used and in order to simulate an infinitely 
long system periodical boundary conditions (Born- Von Karman) were introduced. 



III. GROUND STATE CONFIGURATIONS 



A. Phase Diagram 

The charged particles crystallize in a certain number of chains. Each chain has the same density resulting in a total 
one-dimensional density ng. It is then possible to calculate the energy per particle for each configuration and to check 
the favored one as a function of the parameters of the system. If a is the separation between two adjacent particles 
in the same chain, we can define the dimensionless linear density rie = Ir^/a, where I is the number of chains. 

In the case of multiple chains, in order to have a better packing (or in other words to minimize the interaction 
energy by maximizing the separation among particles in different chains), the chains are staggered with respect to 
each other by a/2 in the x direction. In an infinite lattice this will lead to the hexagonal WC 20j. We calculated the 
energy per particle as a function of the density for the first six possible configurations of the system. 

If the particles crystallize in a single chain, the minimum energy is obtained when the particles are placed on the 
y axis, where the confining potential is zero. In this case the linear density is rie = ro/a and the x coordinate of the 
particles are Xi = ia, with i = 0,±1,±2, ...,±oo. The energy per particle is : 



oo ^ 

El = ne^- exp(-Kj7rie). (3) 



The case of Coulomb interaction is treated using the Ewald summation method so that the summation over long 
distance can be done effectively. Following the standard procedure ,2^] we obtain for Ei : 

Ei{k = 0) = ^?zm[V2e-2-J-$i(jV2He) + V$2(n^(a: - jf) + ^^2{nlx^) - ^-], 



(4) 

where $i(a;) = y/n dt exp[-t^) l/t, <^2{x) = eTk{y/x) and erfc(2/) = I - ^ e~^' dt. 

The first summation contains a divergent term at j = coming from the lower limit of the integration in the 
function $i(a; = 0). This divergence is remedied if we subtract the interaction energy {E},) of the negatively charged 
particles with the positive background which also diverges logarithmically in one dimension. In that case we can 
proceed using the limit limx~^ erf(x) = 2/y/Tr: 
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AEi = Ei{n = 0)~Eb^ ^[V2$i(jV2He) + V$2(H2i2)] _ -l^n^. (5) 



In the two-chain configuration the particles crystallize in two parallel lines separated by a distance d and displaced 
by a distance a/2 along the y-axis. The energy per particle in this case is: 



= 4 + ^y-exp(-2.,7Hg) + ^e^exp(-2.V(,-l/2)^ + cVng) 

where rie = 2ro/a and c = d/a. The first term in (4) is the potential energy due to the confining potential, the 
second term is the energy due to the intra-chain interaction and the last term represents the inter-chain interactions. 
Minimizing E2 with respect to the separation between the chains, c, we obtained the ground state energy for the 
two-chain configuration. 

Similar straightforward but tedious calculations were done for the other multi-chain structures. By symmetry there 
is one intrachain distance in the three-chain structure, two in the four- and five-chain structures and three in the 
six-chain structure. The corresponding expressions for the energy are relegated, for completeness, to Appendix A. 

Calculating the energy minimum for each configuration for different values of at fixed k, we obtain the energy 
per particle E. In Fig. 1 we show E as & function of the density Uf. for k = 1. Notice that for certain density ranges 
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more than one configuration can be stable (this is made more clear in the insets of Fig. 1 for around 2 and 4.7). 
In the low density limit the energy per particle is given by the first term of Eq. (3) : E = ng^exp{—K/nf,), while the 
rest of the curve can be fitted to = -0.0194ng + 0.720ne - 0.245 with an error less than 2.3%. 

Calculating the energy minima for different He and different k we obtain the zero temperature phase diagram of 
Fig. 2. For k = we recover the Coulomb limit. We found that the energy obtained by the analytical method is in 
excellent agreement with the one obtained by our Monte Carlo simulations with a difference between them less than 
0.3 %. 

We observe the following sequence of transitions as the density increases : from one to the two-chain structure 
then to the four-chain configuration, back to the three-chain and again to four and then to five, six etc. Notice the 
remarkable fact that between the 2 and 3 chain configuration there is a small intermediate region where a four-chain 
configuration has a lower energy. For all other transitions the number of chains increases only by one unit, i.e. 
n n + 1. The relative lateral position of the different chains are depicted in Fig. 3 as a function of the density rig. 
In the case of two and three-chain the inter-chain distance increases as the density increases. This is also true for the 
four-chain configuration too, with some differences. In the first four-chain regime of the phase diagram, the distance 
between the two internal chains is larger than the distance between the internal chains and the external ones, in the 
second regime the behavior of the system is the opposite with the distance between internal and external chains larger 
than the one between internal chains. For the other structures the interchain distance is always a growing function 
of the density. It is evident that only the first transition is continuous with a clear bifurcation. 

In order to gain some insight on the distribution of the energy in this anisotropic system we present in Fig. 4 the 
energy per particle for each chain. This is computed by considering a particle at a particular chain and taking into 
account all the interactions with the rest of the particles. The cases of interest are the configurations for which it 
is possible to distinguish internal from external chains and may be related to the difference in the melting behavior 
which is discussed in Sec. V. The interesting observation is that in every case the energy per particle is larger in the 
external chain than the internal ones. 

This asymmetry reflects the fact that for each particle residing in an external chain the gain in energy due to the 
confining potential is higher than the difference in the Coulomb energy due to the lack of symmetric neighboring 
chains, as compared to a particle residing in an internal chain. E.g. for a three-chain system where the middle chain 
is the 0th and the external ones are denoted by +1 and -1, we have for the energy of two particles : 

E±l — Ef) ~ Econf,±l + Ecoulomb, + l-l " Ecoulomb,±l,0 > 0, (7) 

where T!^couiomb,a,(3 denotes the Coulomb energy of a particle residing in chain a interacting with the particles in chain 
(3 and Econ/,a denotes its confining energy. 

In the case of the first density regimes where the four-chain structure is optimal this difference is not large due 
to the fact that the internal distance is less than the external. On the contrary, the difference is much larger in 
the second regime of the four-chain structure. Another interesting observation is that as we approach the limit of 
Coulomb interactions (k ^ 1) the energy difference tends to vanish and the system behaves isotropically. 

B. Structural transitions 

We have seen that by increasing the density, the system changes its configuration, in other words it undergoes a 
"structural transition". It is a natural question to study the order of these transitions. For this purpose the derivative 
of the energy with respect to the density was calculated which is shown in Fig. 5 for the case of k = 1. Only 
the transition between the one and the two-chain configuration is continuous and all the others are discontinuous. 
This conclusion agrees with the results of Fig. 3 where discontinuous changes of the lateral position of the particles 
correspond to first order transitions. The transition 1 ^ 2 is a "zig-zag" transition (Fig. 6). The transition 
2^4 occurs through a "zig-zag" transition of each of the two chains accompanied by a shift of a/4 along the chain, 
which makes it a weakly discontinuous transition (Fig. 6). In principle, these kind of almost "zig-zag" transitions are 
possible for three-, four-, five- and six- chains to result into six-, eight-, ten- and twelve- chain structures, respectively. 
Actually, these were observed during the numerical simulations, especially for very small value of k, but they represent 
metastable states and are not the most energetically favored configurations. 

C. Limit of short range interaction and large density 

In order to make the connection with the regime where the hard core potential can be used as a working hypothesis, 
we investigate the limit A <C a. It can be shown that the variation of the distances between chains can be neglected 
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and in the limit where mujQ^W'^ ^ /{ae) (W is the width of the strip), following the spirit of the hydrodynamic 
consideration of Koulakov and Shklovskii (2^J the difference in the distance between chains at the borders d{±W/2) 
and at the center d{0) follows the relation: 

Sdo = do{W/2) -do{0) ^ Xlnl, (8) 



where I is the number of chains and do — d'^ + a? / A. 

This can be estimated by considering the pressure Uyy in the crystal exercised by the external potential. Adopting 
a method similar to Ref . |25L l26l| : 

(Jyv^~S{a)^^{—-y ), (9) 

where 3/4 < S{a) < 1 and a is the Poisson ratio. We assume a uniform density n and S{a) w 1. Then, balancing the 
force by the pressure and the interaction forces we get (in this estimate we keep the dimensions for clarity): 

'''' -exv{-4^TW^I\)-''-^{W-lA-y\ (10) 



e(d2 + a2/4) 



from this relation 



4q2 



do(y) ^ ^d[yf + aV4 ^ A/n[ ^^^^,^^^,^^ _ ^ ^,^^^ ], (11) 

subtracting the values of do at — a/2 and we obtain Eq. (6). 

Therefore in the case of very short-range interaction bd <^ d(0). Then one can adopt the hard core potential and 
essentially the total energy becomes the sum of the energy of each particle due to the confining potential. The average 
energy per unit length E/L then reads : 

1 

E/L^-mu:^ — . (12) 
24 La 

IV. NORMAL MODES 
A. Normal modes in the absence of an external magnetic field 

We now turn to the calculation of the normal modes of the system, following the standard harmonic approximation 
[SOj l and exploiting the translational invariance of the system along the x direction. The number of chains determines 
the number of particles in each unit cell and therefore the number of degrees of freedom per unit cell. So if / is the 
number of chains there will be 21 branches for the normal mode dispersion curves: / acoustical branches as well as I 
optical ones. Note that for ordinary bidimensional crystals there are 2 acoustical branches and 2r — 2 optical branches, 
if r is the number of atomic species in the unit cell. We present the results for the one, two and three-chain structure 
in Fig. 7. Note that for the one chain structure the unit cell consists of a single particle, i.e. r — 1, and therefore one 
expects only a single acoustical branch and no optical branch. The appearance of the optical branch is a consequence 
of the presence of the confining potential in the y-direction. Note that for fc ^ 0, uJopt ~ which corresponds to 
the center of mass motion of the system in the confining potential. 

In order to find the eigenmodes we solve the system of equations : 

[U}'^5ap,ij — Dal3.ij)Ql3,'j — 0, (13) 

where Qp.j is the displacement of the particle j from its equilibrium position in the (3 direction, (a, (3) = (x, y), 5ap,ij 
is a unit matrix and Dap^ij is the dynamical matrix defined by: 

m ^ — ' 

V 

where v is an integer assigned to each unit cell and the force constants are : 



^ix-x')^ + {y-y')^ 
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evaluated at x — x' G {aiy, a{v + 1/2)}, y — y' = relevant interchain distance. And 

</'a,/3(i' = O) = -^0„,^(z.). (16) 

All the frequencies are measured in unit of ujq/^/2. In Appendix B we present for completeness the expressions for 
the matrix where as an example the modes for the three-chain structure were calculated. 

The main feature is the softening of the optical mode of the one chain structure at the values of ng and k where 
the structural transition is observed ( "zig-zag" transition) accompanied by a hardening of the acoustical branch (Fig. 
8), which confirms that 1 —s- 2 is a continuous transition as asserted before. 

Studying the eigenvectors of the dynamical matrix it is easy to recognize that the optical modes are identified with 
the motion in the direction of confinement (y direction) while the acoustical modes are identified with the motion in 
the unconfined x direction. 

The eigenfrequencies for the single chain are given by : Uac = for the acoustical branch and ujopt = VI + ^2 
for the optical branch, where Ai and A2 are defined in Appendix B. 

In the limit of small wavenumbers k, the summations can be done analytically and we obtain: 

u^oAk) = {1 - [-ln{l - e~-rn.) + ^ Yf^iea'Y/^ (18) 

He 1 — e 

which gives explicitly the dependence of the modes on the density and the screening parameter. In the limit n/rie 3> 1 



u:ac{k) = e"-/"=-^|A:|a, (19) 



u^opt{k) = l-e-^/^'^ea\ (20) 



while in the opposite limit n/ue <C 1 



u^acik) = [l + lni^)]^/'nl/^\k\a, (21) 

2 K 



cooptik) - {l-[l + lni^)]nieaY/-\ (22) 



There is a remarkable difference in the optical branch of the spectrum between the single chain and the two and 
three-chain structures. In the first case the frequency of the optical branch decreases as the wavenumber k increases, 
while for the two and three-chain structures the optical frequency increases. In the single chain configuration the 
optical mode corresponds to oscillations of the particles in the confined direction (see e.g. Fig. 9(b)) which reduces 
the Coulomb repulsive energy. For the two chain configuration the normal modes are shown in Fig. 9(c-g). In fact 
this branch is nothing else than a transverse acoustical mode while the acoustical branch corresponds to longitudinal 
motion Ii.l37ll. 



B. Normal modes in the presence of an external magnetic field 



We now consider the effect of applying a constant magnetic field B in the z direction. For quantum particles, the 
magnetic field can localize the charged particles into cyclotron orbits, therefore aiding the formation of a Wigner 
crystal in the presence of a magnetic field. It is known 271 that in a classical system an external magnetic field 
does not alter the statistical properties of the system and consequently the structural properties and the melting 
temperature are insensitive to the magnetic field strength. But on the other hand the character of motion of the 
particles is altered significantly when the cyclotron frequency is larger than the eigenfrequencies of the system. The 
magneto-phonon spectrum of an infinite 2D Wigner crystal in a magnetic field was obtained in [l^- In the 
presence of B, the system of equations is modified to: 



{uj'^Sap.ij — Dap.ij + i0JLUc£,apSij)Q/3j = 



(23) 
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where ^q./3 is the Levi-Civita tensor and U)c = qB/mc is the cyclotron frequency. In Fig. 10 we show some typical 
dispersion curves for the one and three-chain structures for different values of ujc- It is interesting to notice how the 
optical modes couple with the magnetic field, the optical frequencies follow the cyclotron frequency and for very high 
field strength there is no significant difference between oJopt and ujc- The acoustical frequencies on the other hand 
decrease with the magnetic field strength. For the single chain the eigenfrcqucncies are modified to: 

u;(fc) = {1(1 + + ^2 + c^c') ±\[{l + A,+A.2+ u,^f - 4Ai(l + ^2)]'/'}'/' (24) 

where Ax and A2 are given in Appendix B. For very large field when oj^ 3> {^1,^2,1} the gap between the optical 
branches and the acoustical ones approaches Wc- The optical frequency reflects the cyclotron motion of the system 
which supresses any soft excitation. As before, it is interesting to study the normal modes at the critical density of 
the transition from the one-chain to the two-chain structure (Fig. 11). We observe that there is always softening 
at the same density, independently of the strength of the magnetic field, but with a main difference that for zero 
magnetic field strength the modes which soften is the optical one while when the magnetic field strength is nonzero, 
the acoustic mode is the one that softens. The magnetic field induces a coupling between the acoustical and the 
optical modes and there is an anticrossing between the two branches. Although these findings confirm the previous 
assertion that the presence of B does not alter the structural properties of the system it also reveals the differences 
(softening of the acoustic mode at the same density, influence on the gap between optical and acoustical branches and 
on eigenfrequencies within each branch) which are induced by the magnetic field. 

V. MELTING 

A. General discussion and results 

In this section we study the melting of the WC by MC simulations. After the ground state configuration was 
achieved as explained in Sec. II, the system was heated up by steps of size AT, typically AT = 5 x 10^"*, and 
equilibrated to this new temperature during 10^ 10^ MC steps. In Fig. 12 we show typical trajectories of particles 
as they arise from our MC simulation. It is evident that there is a different behavior of the system in the x and the 
y directions as may be expected by the anisotropy in the two directions. In order to quantify the observations, we 
studied first the potential energy as a function of temperature (Fig. 13). In the crystalline state the potential energy 
of the system increases practically linearly with temperature and then exhibits a very fast increase in a small critical 
temperature range after which it starts to increase linearly again but now with a slightly larger slope. In the latter 
region the system is in the disordered (i.e. liquid) phase. The fast increase of the potential energy is indicative of 
the melting of the WC. To find the critical temperatures we studied, following the spirit of Ref. ||31|, the modified 
Lindemann parameter Lp — {v?)/dr'^, where (w^) is defined by the difference in the mean square displacements of 
neighboring particles from their equilibrium sites and dr is the relevant interparticle distance as we discuss below. 
The quantity (u^) can be written as: 

= ]^(E— E[(^^' - ^^oO - ^^o:'-)]') (25) 

where ( ) means the average over the MC steps, N is the total number of particles in our simulation unit cell and the 
index j denotes the Nnb nearest neighbors of particle i. In order to describe more accurately the difference between 
the two directions, we studied separately (u^) and (u^) as function of temperature. For the melting along the x 
direction, the distance dr is the interparticle distance a introduced in Sec. I while for melting along the y direction 
dr is the interchain distance which is a function of the density rie- 

At low temperatures, the mean square relative displacements slowly increases linearly with temperature as a conse- 
quence of harmonic oscillations of the particles about their equilibrium positions (see Fig. 14). From Fig. 14 we notice 
clearly that this linear increase is larger in the unconfined direction than in the confined direction. In some critical 
temperature region, (u^) and (u^) start to increase very rapidly which is the consequence of the fact that the particles 
have attained sufficient thermal energy that they can jump between different crystallographic positions. According 
to the modified Lindemann criterion MLC, when Lp reaches the (semi-empirical) critical value 0.1 the system melts. 
This criterion was used to define the melting temperature Tm- 

From the corresponding analysis two different melting temperatures, Tx and Ty, can be assigned. The results are 
summarized in the phase diagram of Figs. 15(a-c) for k =0.01, 1 and 3, respectively. There are several interesting 



8 



features in these phase diagrams: (a) the nearly Coulomb system [k = 0.01) has a melting temperature which is on 
average 15-20 % higher than for the screened Coulomb inter-particle interaction with k — I, which has on its turn 
an average melting temperature about 15 % higher than the screened Coulomb system with k ~ 3. Therefore, we 
conclude that the effect of screeening is to reduce the melting temperatures ; (b) a reentrant behavior is observed as 
a function of density, the minima of the melting temperatures occur at the values of the density where the structural 
phase transitions were predicted (see Fig. 2); (c) there is a regime close to each structural transition point where 
the system is frustrated, in the sense that it fluctuates between the two structures. In this regime, which we term as 
frustration regime, the system makes continuous transitions from one metastable state to the other which strongly 
reduces the melting temperature; (d) for n — 1 and k = 3, there is a region in density for which the system melts 
first in the unconfined direction while it is not melted in the confined one. This regime resembles the findings of 
Ref. |l5l | in the regime termed as locked floating solid. For the Coulomb limit there is no evidence of anisotropic 
melting within the error bars of our simulation. The system behaves more isotropic; (e) the first four-chain regime 
(see insets of Figs. 15a-c) is unstable with respect to temperature fluctuations as it is reflected in the relative low 
melting temperature. In this region, melting occurs first in the confined direction as a consequence of the particular 
structural properties -the distance between the two internal chains is larger than the distance between an internal 
chain and the adjacent external one- which makes the system unstable in the y direction. In the rest of the diagram 
there is evidence that the melting either starts from the unconfined direction (e.g. it is clear in the single chain and 
in the low density limit of the two-chains) or the system melts simultaneously in both directions. (/) the single-chain 
structure shows a relatively large melting temperature as obtained by the MLC and deserves more attention. The 
study of the single-chain is therefore postponed to the next subsection. 

Furthermore, note that the MLC only takes into account the displacement of the particles relative to the position 
of their neighbors and consequently is only a measure of the local order of the system. 

Another natural question that arises is whether there is anisotropic melting with respect to external and internal 
chains in the multi-chain structures or in other words if melting starts from the edges as observed in the experiment 
of Ref. 01 with electrons on liquid helium. The number of filaments that were observed in the experiment was 
approximately 20; we have simulated the trajectories of some multi-chain structures and the results are presented 
in Fig. 18. In this picture it is clear that the most external chains are already melted while the internal ones 
are still ordered. Edge melting has also been demonstrated in the presence of a strong magnetic field in Ref. [s^ 
using Hartree-Fock calculations in a two-dimensional Wigner crystal with edges. With the aid of many numerical 
simulations of multi-chain systems at different densities we observed that this kind of melting is present in our system 
when the density is close enough to the critical density of a structural transition. Close to the structural transition 
many metastable states appear with a different number of particles per chain, that is in the most external chains there 
are less particles than in the internal ones. Thus the particles at the most external chains have larger displacements 
from their equilibrium positions in order to attain the stability of the structure. Furthermore, we calculated the 
average root mean square displacements of the particles from their equilibrium position chain by chain and also (w^) 
chain by chain and we actually noted that these quantities are slightly larger for external chains at temperatures 

below the critical one. In Fig. 19 we present the temperature dependence of the standard deviation Sx'^ = 
{{ux — (ux))^) and Sx"^ = {{uy — {uy)Y) for the external and internal chains in the four-chain structure. It is evident 
that the position of the particles at the edges fluctuates substantially more than the particles at the interior. We can 
conjecture that, according to this physical picture, melting can start from the edges. However, for up to the six-chain 
configuration for each chain the quantities (m^) reached the critical value, approximately, all at the same temperature. 
Probably, going to a larger number of filaments one can well appreciate a different melting temperature for external 
and internal chains. Finally, the chain configuration as well as the melting which starts from the direction of the 
chains is supported also by molecular dynamics simulations of the flow of electrons in QID channels [3^. 



B. Melting of the single-chain 

In Figs. 15(a-c) we observe a rather high melting temperature in case of the one-chain structure. The origin of this 
behavior can be traced back to the fact that the MLC takes into account a larger contribution from jumps of particles 
between crystallographic positions which for the single chain structure occurs only at extremely high temperature. 
For the single-chain case the jumps can only occur along the chain which requires a larger energy than jumps of 
particles between different chains. 

To have a better insight we investigated the behavior of Lp for different densities (Fig. 18). We notice that in the 
low density limit (see Fig. 18(a)), Lp w 0.1 is reached in a region in which there is only a gradual increase in (w^) 
which is very different from the multi-chain case (see Fig. 14). Furthermore, (it^) exhibits a sublinear temperature 
increase. 

This calls for the use of other possible criteria in order to clarify the situation. On the other hand if the density is 
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relatively high (see Fig. 18), a fast increase is observed signaling a clear melting of the system. The transition from 
a low temperature linear to sublinear behaviour occurs for Uf, « 0.4. 

To shine light into the posed questions we studied also the pair correlation function at different densities and 
temperatures, as defined by 

9i^) = ^J2<^i^-^^^-''i'^^>' (26) 

where in the summation over N particles in a system of length L, the diagonal terms {i ~ j) are excluded . The 
results are reported in Fig. 19. It is rather evident that the melting temperature is substantially smaller than the 
one obtained from the MLC. In order to better quantify the melting temperature for the one-chain structure we 
investigated the height of the first and second peak of the pair correlation function as function of temperature (see 
Fig. 20) in order to look for a structure or an anomalous jump (as found in Ref. [s^) that could identify the critical 
temperature. As is apparent from Fig. 20 we do not find any abrupt changes. The first and second peak as a function 
of temperature can be fitted by : = a{T/TQ)^^ , where z = {1 or 2} denotes the peak (see the curves in Fig. 20). 

The values of (a,/3) are (2.922, 0.274) for i = 1 and (1.895, 0.216) for i = 2 when = 0.2, (9.320, 0.433) for i 1 
and (7.345, 0.466) for i = 2 when He = 0.5, (14.788, 0.473) for i = 1 and (11.552, 0.501) for z = 2 when Ue = 0.8 and 
in each case the error is less than 1%. 

From the study of the pair correlation function we conclude that at moderate (n^ < 0.2 for k = I) densities, the 
chain is melted at arbitrarily weak temperature. For higher densities the chain retains correlations up to higher values 
of the temperature but these values are less than those obtained by the MLC. 

We noticed from the high density regime (ne > 1), where we reach the multichain structure, that another semi- 
empirical criterion can be formulated using the pair correlation function. If we consider the ratio of the height of 
the fifth peak in g{x) above 1 {H^ — 1) to the height of the first peak above 1 {Hi — 1) at those densities where the 
maximum melting temperature is obtained for the two, three and four chain structures, melting occurs when : 

^^(r™)«0.15 (27) 

Employing this criterion (termed as pair correlation function criterion or PCFC) we obtain the results of Fig. 21, 
were we present both the relevant temperatures obtained by MLC and PCFC. 

It is worth noticing that this criterion does not work well at temperatures close to the structural transitions. The 
reason is that although particles "jump" to new sites in order to attain the new positions, the pair correlation function 
still measures correlations at certain distances and, most importantly, the height of the first peak is substantially 
reduced which artificially enhances the ratio (25). 

Thus the value of 15% which works far from the structural transitions is too high for the regime close to the 
structural transitions. It is therefore evident that the two criteria can work in a complementary manner. 



VI. DISCUSSION AND CONCLUSIONS 



The structural phase transitions and the melting can be studied experimentally using parabolically confined col- 
loidal particles or dusty plasmas in the case of a screened Coulomb inter-particle interaction. Another important 
experimental system are electrons floating on liquid helium where it is possible to achieve relatively narrow QID 
channels on very stable suspended helium films over structured substrates [3^ . Assuming a semicircular profile of 
the liquid surface across the channels then the confining potential is parabolic near the bottom with ujq = {eE'^ /mR) 
[33 where E*^ is the effective holding electric field in the case of the substrate and R the radius of the semicircular 
profile. Assuming a radius of approximately 5 /im, a typical value for E*j^ « lOkV/cm then ujo « 10"'^^ Hz. This in 
turn produces a Tq ~ 60^^. The melting temperatures which have been obtained in the present work are of order 
10^^ X To which results in a melting temperature « 0.5 — IK, a temperature range which is routinely achieved in such 
experiments. Assuming an interelectron distance of approximately 0.1-1 fim leads to a dimensionless linear density 
fig « (0.5 — 3)1, where I is the number of chains. The dilute limit gives the same fie as the one investigated in the 
present work. 

Another issue connected with melting, which deserves interest, is the appearance of topological defects so that a 
KTHNY ^3 scenario of melting is possible. In Ref . [25ll3lll3^ this question was considered in the case of a circular 
confining potential with a finite number of particles. In the case of short range interactions the defects are pushed 
to the surface due to the large price for elastic deformations while in the Coulomb case shear and Young moduli are 
relatively small [25|. 
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Moreover, because of the incommensurability of the circle with the hexagonal Wigner crystal the defects do not 
reside exactly at the borders but in a zone few lattice spacings inside the crystal. Therefore three different regimes 
with different melting temperatures can be detected |38i|. In our case there is no such incommensurability and the 
edges can accommodate the defects easily. This has also been discussed in the case of a quantum Hall bar by Nazarov 
in Ref.lH. 

In conclusion, wc investigated the structural, dynamical properties and melting of a classical quasi- ID system of 
particles interacting through a Yukawa-type potential in the range from Coulomb to very short range interaction in 
the case where the confinement is modelled by an external parabolic potential. The structural transitions are of first 
(primarily) and second order. The normal modes of the system were calculated in the presence of a perpendicular 
magnetic field. In certain regions of the parameter space, there is evidence that melting starts first in the unconfined 
direction before the system melts along the chain direction. Furthermore, we found that shows a reentrant behavior 
as a function of the density of the system and a regime of frustration around each point of structural transition can 
be identified. In the case of the single chain structure, we device a new criterion in order to take into account the 
correlations at different temperatures. The present study is suitable to describe colloidal particles, dusty plasmas and 
electrons floating on liquid helium. 
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VII. APPENDIX A 



The expressions for the energy in the configurations beyond the two-chain structure are presented below. All the 
distances are in units of the interchain distance a between adjacent particles. 
For the three-chain structure : 

/ o /- ^ / exp{~3Ky/ [m ~ 1/2)^ + Ct? jue) 

^3 = exp{~3mK/nej/m + — ^ 



2 



2ne ^ exp{~3K^/m'^ + Acs'^/ue) 03^ neexp{-6c3K/ne) 

' 9 ,h v/-^ + 4c32) + 18c3 ' ^''^ 



where the intrachain distance C3 is a variational parameter. 
For the four-chain structure : 



exp{—4:Km/ne) exp(— 4k-\/ (m — 1/2)^ + (c4 — fiY /Uf^) 



E — exp{-4:Km/ne) _|_ ^ 

A 2-^ rn. 4 ^ 



ne exp(-4K^(m2 + (04 + U? jn^) exp{~An^ (to - 1/2)2 ^_ 4/^2/^^) 



Vm2 + (c4 + /4)2 8 ^(m- 1/2)2 + 4/4' 



Tie exp{-AK^J (to - 1/2)2 ^ 4c42/ng) 8C42 „ J, 2 /~2 

, rie ea;p[-4(c4 fA)n/ne\ , . 

8 C4 -I- 74 

where C4 is the distance of an inner chain and is the distance of an outer chain from the middle of the structure. 
These distances are two variational parameters which have to be optimized numerically. 
For the five-chain structure : 
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10c5^ , 10/5^ Tie exp{-5Km/ne) , 4ne exp{-5K^ (m - 1/2)^ + (cs - /5)^/ne) 

-^5 = ^?3- + + — ^ 



25 £ 1/2)2 + (c5-/5)^ 



4ne e.Tp(— 5K\/mM~C5^/ne) 2ne exp{—5c5K/ne) 4ne exp{—hK\/ {m — 1/2)^ + (05 + f^Y /ue) 



25 



25 



C5 



25 



m=l 



V(m - 1/2)2 + + 



2ne exp{-bK^/m:^ + 4c5Vne) ne exp{-lQc5K/ne) An^ >^ exp{-5K^ (m - 1/2)2 + /g^/^^) 



C5 



m=l 



^(m- 1/2)2 + /5 



2ne ^ e.Tp(-5KY'm2 + 4/5 )/ne) neea;p(-10/5K/ne) 



25 



Vto2 + 4/5 



50/5 



(30) 



with the variational parameters C5 and /s, which are the distance of an inner chain and outer chain respectively from 
the middle of the structure. 
For the six-chain structure : 



12/6^ ^ 12^(62 ^ 12/i6^ I !!£ exp{-6K-m/6) Jh^ , 



nl nl nl '6 



18 



exp{-6K^ (m - 1/2)2 ^ 4/ig2/^^) 
^(m- 1/2)2 + 4/16^) 



+ 



18 



expi-6K^ (m - 1/2)2 + 4/6 Vrte) ea;p(-6Kv/(m - 1/2)2 + 4g6Vne) 



m=l 

00 



(m - 1/2)2 + 4/jg2-) 



m=l 



Tie ■y ea:p(-6KA/ (m - 1/2)2 _|_ (/^^ _ g^/ue) ne 



m=l 



V(m- 1/2)2 +456^ 

exp(-6K V(m - 1/2)2 ^ (/.^ /n^) 
V(m - 1/2)2 + (/^g + /g)2) 



ne -y exp{-6K^{m - 1/2)2 ^ _ /6)2"/ne) ?ie e^p(-6K;-y/m2 + (/l6 - UY)/ne) 
« ^ 1/2)2 + (56-/6)2) q 



Vm2 + (/i6 - /6)2 



_l_ Tie exp{-&\hQ - felK/ue) ^ ^ ^ ea;p(-6KA/ ni^ + {lit; + gQ)'^/ne) ^ Ue exp[-6{h6 + g6)i^/ne 



18 

~ 00 



\h6-f6\ 



^/m^ + {he + gey 



18 



{he + ge) 



ne exp{-6Ky/m'^ + (gg + Uy/rie) ^ Ue exp[-6{ge + f(i)K/ne 



18 



(S6 + /e) 



(31) 



with 96 and /ir the three chain distances from the middle of the crystal, starting from the inner one which are the 
variational parameters. 

VIII. APPENDIX B 

We present the matrix uP'l — D where I is the unit matrix and D is the dynamical matrix, for the calculation of 
the normal modes for the three-chain structure. It reads ; 



a;^I-D = 



/a;2 


-^1 





-As 





-As 





\ 


(^2 


-wo^) 


-A2 


-Ai 





-^6 









a;2-Ai 





-As 












-Ai 





(a;2 - wo^) - A2 











-A 





-As 





a;2_Ai 







V 










-A4 





(a;2-a;o2) 


-A2) 



(32) 



where the parameters are 
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OO ^ 

^1 = fil^^;^exp(-3jK/ne)[2 + 6jK/ne + 9fK'^/nl]{l-cos(kTrj)), 



OO ^ 

'^eX! 97^^^^(~^-^'^/"^)(^ + 3jK/ne)(l - cos{knj)), 

3 = 1 ^'^ 



^3 = "eE o^ ./o^2 e^P(-3^V(j + l/2)^ + cV»e)[co.(fc^(j + l/2))-l](j + l/2)^x 



9^V(j + l/2)^ + C3^ ^ 9^^[0- + l/2)^ + C3^] ^ 3j _ ^ 3«y(jTW + ^^j^^. ^ ^^2)2 + ^3 
ne rig He 

00 ^ 

M = ^eE o^ 1/o^9 exp(-3A.V(j + l/2)^ + C3V^e)[cog(fc^(j + l/2))-l]c3^ X 

,=1 27VU + 1/2)- + C3^ 



+ 1/2)2 + C32 + + + + 3] - (1 + W(Z±W±£E)[(, + 1/2)2 ^ 

00 ^ 

^5 = - exp{~?>K\/ p + 4c2/ne)[cos(fc7rj) — 1] x 

^6 = "-eY] y } exp{-3K^/f + 4c32/ne)[cos(fc7rj) - 1] X 

^ 27v/j2 + 4C32 



2 



Tie He 

where k = ak/n, is the dimensionless wavenumber. 

The modes for the single-chain are obtained by the top left part of the matrix (32) which forms a 2 x 2 submatrix 
and involves the elements Ai and A2 which have exactly the same form with the substitution fie/S — > fie- 

Similarly, the modes for the two-chain structure can be obtained by the 3x3 submatrix which is included in the 
top left part of matrix (32) and involves the elements Ai, A2 and A3 with the substitution he/3 h^jl. 
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FIG. 1: The energy per particle as a function of (iensity for k=1. 




FIG. 2: The zero temperature structural phase diagram. 




FIG. 3: The lateral position of the chains in the Wigner crystal state as a function of the linear density for n = 1. 
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FIG. 5: The derivative of the energy with respect to the density for a = 1. Only the transition from one to two-wires is 
continuous (second order) the rest are first order. 
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FIG. 6: The mechanism of the structural transitions 1 — > 2 chains and 2 — * 4 chains. 
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FIG. 11: The magnetic field dependence of the softening of the phonon mode at the structural transition from one to two-chains. 
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FIG. 12: Particle trajectories for 3 x lO'^ MC steps for three different temperatures and three different values of the density at 

K = l. 
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FIG. 14: The mean square relative displacements for k= 1, ne= 1.6 in the case of the two-chain configuration. The dashed 
horizontal line corresponds to the modified Lindemann criterion (MLC) in the unconfined direction, while the dotted line 
corresponds to the MLC in the confined direction. 
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FIG. 17: Temperature dependence of the displacements in both the unconfined and confined directions for external and internal 
chains in the four-chain structure. 
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FIG. 18: The behavior of the Lindemann parameter, for the single-chain regime at four different densities. It shows how the 
hnear regime at higher densities becomes subUnear at lower densities. 
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FIG. 19: The pair correlation lunction at different temperatures, for three different densities, for the single chain configuration. 
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FIG. 20: The height of the first and second peaks of tlie pair correlation function for the single chain as a function of temperature 
for three different densities. The hnes axe the best fits with the function a{T /To)~^ . 
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